﻿Imports System.Drawing

Namespace Maps.Conv
    'Public Class PointF
    '    Public x As Double
    '    Public y As Double
    '    Public Sub New(ByVal x As Double, ByVal y As Double)
    '        Me.x = x
    '        Me.y = y
    '    End Sub

    'End Class
    Public Class WGStoUTM
        Const UTMScaleFactor = 0.99960000000000004
        Const sm_a = 6378137.0
        Const sm_b = 6356752.3140000002
        'Public Shared Function ToUTM(ByVal p As PointF) As PointF
        '    'Dim xy(1) As Double

        '    'Dim lng = p.X : Dim lat = p.Y
        '    Dim lat = p.Y : Dim lng = p.X

        '    ' Compute the UTM zone.
        '    Dim zone As Integer = CInt(Math.Floor((lng + 180.0) / 6) + 1)
        '    zone = 32

        '    Dim xy = LatLonToUTMXY(Utils.Math.DegToRad(lat), Utils.Math.DegToRad(lng), zone)

        '    Return New PointF(xy(0), xy(1))

        'End Function
        Public Shared Function ToUTM(ByVal p As PointF, Zone As Integer) As PointF
            'Dim xy(1) As Double

            'Dim lon = p.Y : Dim lat = p.X
            Dim lat = p.Y : Dim lng = p.X
            ' Compute the UTM zone.

            Dim xy = LatLonToUTMXY(Utils.Math.DegToRad(lat), Utils.Math.DegToRad(lng), Zone)

            Return New PointF(xy(0), xy(1))

        End Function
        Public Shared Function ToUTMAuto(ByVal p As PointF) As PointF
            'Dim xy(1) As Double

            Dim lon = p.Y : Dim lat = p.X

            ' Compute the UTM zone.
            Dim zone As Integer = CInt(Math.Floor((lon + 180.0) / 6) + 1)

            Dim xy = LatLonToUTMXY(Utils.Math.DegToRad(lat), Utils.Math.DegToRad(lon), zone)

            Return New PointF(xy(0), xy(1))

        End Function

        Private Shared Function LatLonToUTMXY(ByRef lat As Double, ByRef lng As Double, ByRef zone As Integer) As Double()
            'Dim xy(1) As Double

            Dim xy = MapLatLonToXY(lat, lng, UTMCentralMeridian(zone))

            ' Adjust easting and northing for UTM system. */
            xy(0) = xy(0) * UTMScaleFactor + 500000.0
            xy(1) = xy(1) * UTMScaleFactor
            If (xy(1) < 0.0) Then
                xy(1) = xy(1) + 10000000.0
            End If

            Return xy

        End Function
        Private Shared Function UTMCentralMeridian(ByRef zone As Integer) As Double

            Dim cmeridian As Double = Utils.Math.DegToRad(-183.0 + (zone * 6.0))

            Return cmeridian

        End Function
        Private Shared Function MapLatLonToXY(ByRef phi As Double, ByRef lambda As Double, ByRef lambda0 As Double) As Double()
            Dim xy(1) As Double
            Dim N, nu2, ep2, t, t2, l As Double
            Dim l3coef, l4coef, l5coef, l6coef, l7coef, l8coef As Double
            Dim tmp As Double

            '/* Precalculate ep2 */
            ep2 = (Math.Pow(sm_a, 2.0) - Math.Pow(sm_b, 2.0)) / Math.Pow(sm_b, 2.0)

            '/* Precalculate nu2 */
            nu2 = ep2 * Math.Pow(Math.Cos(phi), 2.0)

            '/* Precalculate N */
            N = Math.Pow(sm_a, 2.0) / (sm_b * Math.Sqrt(1 + nu2))

            '/* Precalculate t */
            t = Math.Tan(phi)
            t2 = t * t
            tmp = (t2 * t2 * t2) - Math.Pow(t, 6.0)

            '/* Precalculate l */
            l = lambda - lambda0

            '/* Precalculate coefficients for l**n in the equations below
            '   so a normal human being can read the expressions for easting
            '   and northing
            '   -- l**1 and l**2 have coefficients of 1.0 */
            l3coef = 1.0 - t2 + nu2

            l4coef = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2)

            l5coef = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2 - 58.0 * t2 * nu2

            l6coef = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2 - 330.0 * t2 * nu2

            l7coef = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2)

            l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2)

            '/* Calculate easting (x) */
            xy(0) = N * Math.Cos(phi) * l _
                + (N / 6.0 * Math.Pow(Math.Cos(phi), 3.0) * l3coef * Math.Pow(l, 3.0)) _
                + (N / 120.0 * Math.Pow(Math.Cos(phi), 5.0) * l5coef * Math.Pow(l, 5.0)) _
                + (N / 5040.0 * Math.Pow(Math.Cos(phi), 7.0) * l7coef * Math.Pow(l, 7.0))

            '/* Calculate northing (y) */
            xy(1) = ArcLengthOfMeridian(phi) _
                + (t / 2.0 * N * Math.Pow(Math.Cos(phi), 2.0) * Math.Pow(l, 2.0)) _
                + (t / 24.0 * N * Math.Pow(Math.Cos(phi), 4.0) * l4coef * Math.Pow(l, 4.0)) _
                + (t / 720.0 * N * Math.Pow(Math.Cos(phi), 6.0) * l6coef * Math.Pow(l, 6.0)) _
                + (t / 40320.0 * N * Math.Pow(Math.Cos(phi), 8.0) * l8coef * Math.Pow(l, 8.0))

            Return xy

        End Function
        Private Shared Function ArcLengthOfMeridian(ByRef phi As Double)
            Dim alpha, beta, gamma, delta, epsilon, n As Double
            Dim result As Double

            '/* Precalculate n */
            n = (sm_a - sm_b) / (sm_a + sm_b)

            '/* Precalculate alpha */
            alpha = ((sm_a + sm_b) / 2.0) _
               * (1.0 + (Math.Pow(n, 2.0) / 4.0) + (Math.Pow(n, 4.0) / 64.0))

            '/* Precalculate beta */
            beta = (-3.0 * n / 2.0) + (9.0 * Math.Pow(n, 3.0) / 16.0) _
               + (-3.0 * Math.Pow(n, 5.0) / 32.0)

            '/* Precalculate gamma */
            gamma = (15.0 * Math.Pow(n, 2.0) / 16.0) _
                + (-15.0 * Math.Pow(n, 4.0) / 32.0)

            '/* Precalculate delta */
            delta = (-35.0 * Math.Pow(n, 3.0) / 48.0) _
                + (105.0 * Math.Pow(n, 5.0) / 256.0)

            '/* Precalculate epsilon */
            epsilon = (315.0 * Math.Pow(n, 4.0) / 512.0)

            '/* Now calculate the sum of the series and return */
            result = alpha _
            * (phi + (beta * Math.Sin(2.0 * phi)) _
                + (gamma * Math.Sin(4.0 * phi)) _
                + (delta * Math.Sin(6.0 * phi)) _
                + (epsilon * Math.Sin(8.0 * phi)))

            Return result

        End Function
    End Class


    Public Class UTMtoWSG
        'Public Shared Function ToWGS(ByVal p As PointF) As PointF

        '    Dim f = p.X : Dim f1 = p.Y

        '    Dim j As Double = 32

        '    Dim d As Double = 0.99960000000000004
        '    Dim d1 As Double = 6378137
        '    Dim d2 As Double = 0.0066943799999999998

        '    Dim d4 As Double = (1 - Math.Sqrt(1 - d2)) / (1 + Math.Sqrt(1 - d2))
        '    Dim d15 As Double = f1 - 500000
        '    Dim d16 As Double = f
        '    Dim d11 As Double = ((j - 1) * 6 - 180) + 3

        '    Dim d3 As Double = d2 / (1 - d2)
        '    Dim d10 As Double = d16 / d
        '    Dim d12 As Double = d10 / (d1 * (1 - d2 / 4 - (3 * d2 * d2) / 64 - (5 * Math.Pow(d2, 3)) / 256))
        '    Dim d14 As Double = d12 + ((3 * d4) / 2 - (27 * Math.Pow(d4, 3)) / 32) * Math.Sin(2 * d12) + ((21 * d4 * d4) / 16 - (55 * Math.Pow(d4, 4)) / 32) * Math.Sin(4 * d12) + ((151 * Math.Pow(d4, 3)) / 96) * Math.Sin(6 * d12)
        '    Dim d13 As Double = d14 * 180 / Math.PI
        '    Dim d5 As Double = d1 / Math.Sqrt(1 - d2 * Math.Sin(d14) * Math.Sin(d14))
        '    Dim d6 As Double = Math.Tan(d14) * Math.Tan(d14)
        '    Dim d7 As Double = d3 * Math.Cos(d14) * Math.Cos(d14)
        '    Dim d8 As Double = (d1 * (1 - d2)) / Math.Pow(1 - d2 * Math.Sin(d14) * Math.Sin(d14), 1.5)

        '    Dim d9 As Double = d15 / (d5 * d)
        '    Dim d17 As Double = d14 - ((d5 * Math.Tan(d14)) / d8) * (((d9 * d9) / 2 - (((5 + 3 * d6 + 10 * d7) - 4 * d7 * d7 - 9 * d3) * Math.Pow(d9, 4)) / 24) + (((61 + 90 * d6 + 298 * d7 + 45 * d6 * d6) - 252 * d3 - 3 * d7 * d7) * Math.Pow(d9, 6)) / 720)
        '    d17 = d17 * 180 / Math.PI
        '    Dim d18 As Double = ((d9 - ((1 + 2 * d6 + d7) * Math.Pow(d9, 3)) / 6) + (((((5 - 2 * d7) + 28 * d6) - 3 * d7 * d7) + 8 * d3 + 24 * d6 * d6) * Math.Pow(d9, 5)) / 120) / Math.Cos(d14)
        '    d18 = d11 + d18 * 180 / Math.PI

        '    'Dim xy(1) As Double

        '    'xy(0) = d18
        '    'xy(1) = d17

        '    Return New PointF(d18, d17)

        'End Function

        'Const sm_a As Double = 6378137.0
        'Const sm_b As Double = 6356752.3140000002
        'Const sm_EccSquared As Double = 0.00669437999013
        'Const UTMScaleFactor As Double = 0.99960000000000004

        'Public Shared Function MapXYToLatLon(ByVal x As Double, ByVal y As Double, ByVal lambda0 As Double) As PointF



        '    ' Get the value of phif, the footpoint latitude. 



        '    Dim phif As Double = FootpointLatitude(y)

        '    ' Precalculate ep2 

        '    Dim ep2 As Double = (Math.Pow(sm_a, 2.0) - Math.Pow(sm_b, 2.0)) / Math.Pow(sm_b, 2.0)

        '    ' Precalculate cos (phif) 

        '    Dim cf = Math.Cos(phif)

        '    ' Precalculate nuf2 

        '    Dim nuf2 = ep2 * Math.Pow(cf, 2.0)

        '    ' Precalculate Nf and initialize Nfpow 

        '    Dim Nf = Math.Pow(sm_a, 2.0) / (sm_b * Math.Sqrt(1 + nuf2))
        '    Dim Nfpow = Nf

        '    ' Precalculate tf 

        '    Dim tf = Math.Tan(phif)
        '    Dim tf2 = tf * tf
        '    Dim tf4 = tf2 * tf2

        '    ' Precalculate fractional coefficients for x**n in the equations 
        '    '            below to simplify the expressions for latitude and longitude. 

        '    Dim x1frac = 1.0 / (Nfpow * cf)

        '    Nfpow *= Nf
        '    ' now equals Nf**2) 
        '    Dim x2frac = tf / (2.0 * Nfpow)

        '    Nfpow *= Nf
        '    ' now equals Nf**3) 
        '    Dim x3frac = 1.0 / (6.0 * Nfpow * cf)

        '    Nfpow *= Nf
        '    ' now equals Nf**4) 
        '    Dim x4frac = tf / (24.0 * Nfpow)

        '    Nfpow *= Nf
        '    ' now equals Nf**5) 
        '    Dim x5frac = 1.0 / (120.0 * Nfpow * cf)

        '    Nfpow *= Nf
        '    ' now equals Nf**6) 
        '    Dim x6frac = tf / (720.0 * Nfpow)

        '    Nfpow *= Nf
        '    ' now equals Nf**7) 
        '    Dim x7frac = 1.0 / (5040.0 * Nfpow * cf)

        '    Nfpow *= Nf
        '    ' now equals Nf**8) 
        '    Dim x8frac = tf / (40320.0 * Nfpow)

        '    ' Precalculate polynomial coefficients for x**n. 
        '    '            -- x**1 does not have a polynomial coefficient. 

        '    Dim x2poly = -1.0 - nuf2

        '    Dim x3poly = -1.0 - 2 * tf2 - nuf2

        '    Dim x4poly = 5.0 + 3.0 * tf2 + 6.0 * nuf2 - 6.0 * tf2 * nuf2 - 3.0 * (nuf2 * nuf2) - 9.0 * tf2 * (nuf2 * nuf2)

        '    Dim x5poly = 5.0 + 28.0 * tf2 + 24.0 * tf4 + 6.0 * nuf2 + 8.0 * tf2 * nuf2

        '    Dim x6poly = -61.0 - 90.0 * tf2 - 45.0 * tf4 - 107.0 * nuf2 + 162.0 * tf2 * nuf2

        '    Dim x7poly = -61.0 - 662.0 * tf2 - 1320.0 * tf4 - 720.0 * (tf4 * tf2)

        '    Dim x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2)

        '    Dim philambda = New PointF()
        '    ' Calculate latitude 

        '    philambda.X = phif + x2frac * x2poly * (x * x) + x4frac * x4poly * Math.Pow(x, 4.0) + x6frac * x6poly * Math.Pow(x, 6.0) + x8frac * x8poly * Math.Pow(x, 8.0)

        '    ' Calculate longitude 

        '    philambda.Y = lambda0 + x1frac * x + x3frac * x3poly * Math.Pow(x, 3.0) + x5frac * x5poly * Math.Pow(x, 5.0) + x7frac * x7poly * Math.Pow(x, 7.0)

        '    Return philambda
        'End Function
        'Public Shared Function FootpointLatitude(ByVal y As Double) As Double
        '    ' Precalculate n (Eq. 10.18) 

        '    Dim n = (sm_a - sm_b) / (sm_a + sm_b)

        '    ' Precalculate alpha_ (Eq. 10.22) 

        '    ' (Same as alpha in Eq. 10.17) 

        '    Dim alpha_ = ((sm_a + sm_b) / 2.0) * (1 + (Math.Pow(n, 2.0) / 4) + (Math.Pow(n, 4.0) / 64))

        '    ' Precalculate y_ (Eq. 10.23) 

        '    Dim y_ = y / alpha_

        '    ' Precalculate beta_ (Eq. 10.22) 

        '    Dim beta_ = (3.0 * n / 2.0) + (-27.0 * Math.Pow(n, 3.0) / 32.0) + (269.0 * Math.Pow(n, 5.0) / 512.0)

        '    ' Precalculate gamma_ (Eq. 10.22) 

        '    Dim gamma_ = (21.0 * Math.Pow(n, 2.0) / 16.0) + (-55.0 * Math.Pow(n, 4.0) / 32.0)

        '    ' Precalculate delta_ (Eq. 10.22) 

        '    Dim delta_ = (151.0 * Math.Pow(n, 3.0) / 96.0) + (-417.0 * Math.Pow(n, 5.0) / 128.0)

        '    ' Precalculate epsilon_ (Eq. 10.22) 

        '    Dim epsilon_ = (1097.0 * Math.Pow(n, 4.0) / 512.0)

        '    ' Now calculate the sum of the series (Eq. 10.21) 

        '    Return y_ + (beta_ * Math.Sin(2.0 * y_)) + (gamma_ * Math.Sin(4.0 * y_)) + (delta_ * Math.Sin(6.0 * y_)) + (epsilon_ * Math.Sin(8.0 * y_))
        'End Function

        Public Shared Function ToWGS(ByVal p As PointF, Zone As Integer) As PointF
            'Dim Zone = 32
            Dim k0 = 0.9996D
            Dim a = 6378137D
            Dim f = 1 / 298.2572236D
            Dim b = a * (1 - f)
            Dim drad = Math.PI / 180

            Dim e = Math.Sqrt(1 - (b / a) * (b / a))

            Dim e0 = e / Math.Sqrt(1 - e * e)

            Dim esq = (1 - (b / a) * (b / a))

            Dim e0sq = e * e / (1 - e * e)
            Dim x = p.X

            If (x < 160000 Or x > 840000) Then
                'alert("Outside permissible range of easting values \n Results may be unreliable \n Use with caution");} 
            End If

            Dim y = p.Y
            If (y < 0) Then
                'alert("Negative values not allowed \n Results may be unreliable \n Use with caution");}
            End If


            If (y > 10000000) Then
                'alert("Northing may not exceed 10,000,000 \n Results may be unreliable \n Use with caution");}
            End If


            Dim zcm = 3 + 6 * (Zone - 1) - 180

            Dim e1 = (1 - Math.Sqrt(1 - e * e)) / (1 + Math.Sqrt(1 - e * e))

            Dim M0 = 0

            Dim M = M0 + y / k0

            'if (document.getElementById("SHemBox").checked === true){M=M0+(y-10000000)/k;}

            Dim mu = M / (a * (1 - esq * (1 / 4 + esq * (3 / 64 + 5 * esq / 256))))

            Dim phi1 = mu + e1 * (3 / 2 - 27 * e1 * e1 / 32) * Math.Sin(2 * mu) + e1 * e1 * (21 / 16 - 55 * e1 * e1 / 32) * Math.Sin(4 * mu)

            phi1 = phi1 + e1 * e1 * e1 * (Math.Sin(6 * mu) * 151 / 96 + e1 * Math.Sin(8 * mu) * 1097 / 512)

            Dim C1 = e0sq * Math.Pow(Math.Cos(phi1), 2)

            Dim T1 = Math.Pow(Math.Tan(phi1), 2)

            Dim N1 = a / Math.Sqrt(1 - Math.Pow(e * Math.Sin(phi1), 2))

            Dim R1 = N1 * (1 - e * e) / (1 - Math.Pow(e * Math.Sin(phi1), 2))

            Dim D = (x - 500000) / (N1 * k0)

            Dim phi = (D * D) * (1 / 2 - D * D * (5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * e0sq) / 24)

            phi = phi + Math.Pow(D, 6) * (61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * e0sq - 3 * C1 * C1) / 720

            phi = phi1 - (N1 * Math.Tan(phi1) / R1) * phi



            '//Output Latitude

            Dim lat = Math.Floor(1000000 * phi / drad) / 1000000



            '//Longitude

            Dim lng = D * (1 + D * D * ((-1 - 2 * T1 - C1) / 6 + D * D * (5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * e0sq + 24 * T1 * T1) / 120)) / Math.Cos(phi1)

            Dim lngd = zcm + lng / drad



            '//Output Longitude

            lng = Math.Floor(1000000 * lngd) / 1000000

            Return New PointF(lng, lat)

        End Function
    End Class


    Public Class WGSToPixel

        Private PixelTileSize As Double = 256.0
        Private DegreesToRadiansRatio As Double = 180.0 / Math.PI
        Private RadiansToDegreesRatio As Double = Math.PI / 180.0
        Private PixelGlobeCenter As PointF
        Private XPixelsToDegreesRatio As Double
        Private YPixelsToRadiansRatio As Double

        Public Sub New(ByVal zoomLevel As Double)
            Dim pixelGlobeSize = Me.PixelTileSize * Math.Pow(2.0, zoomLevel)
            Me.XPixelsToDegreesRatio = pixelGlobeSize / 360.0
            Me.YPixelsToRadiansRatio = pixelGlobeSize / (2.0 * Math.PI)
            Dim halfPixelGlobeSize = Convert.ToSingle(pixelGlobeSize / 2.0)
            Me.PixelGlobeCenter = New PointF(halfPixelGlobeSize, halfPixelGlobeSize)
        End Sub

        Public Function FromCoordinatesToPixel(ByVal coordinates As PointF) As PointF
            Dim x = Math.Round(Me.PixelGlobeCenter.X + (coordinates.X * Me.XPixelsToDegreesRatio))
            Dim f = Math.Min(Math.Max(Math.Sin(coordinates.Y * RadiansToDegreesRatio), -0.9999D), 0.9999D)
            Dim y = Math.Round(Me.PixelGlobeCenter.Y + 0.5 * Math.Log((1.0 + f) / (1.0 - f)) * -Me.YPixelsToRadiansRatio)
            Return New PointF(Convert.ToSingle(x), Convert.ToSingle(y))
        End Function
    End Class


End Namespace
